Geminal wavefunctions with Jastrow correlation: a first application to atoms 
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We introduce a simple generalization of the well known geminal wavefunction already applied in 
Quantum Chemistry to atoms and small molecules. The main feature of the proposed wavefunction 
is the presence of the antisymmetric geminal part together with a Jastrow factor. Both the geminal 
and the Jastrow play a crucial role in determining the remarkable accuracy of the many-body state: 
the former permits the correct treatment of the nondynamic correlation effects, the latter allows the 
wavefunction to fulfill the cusp conditions and makes the geminal expansion rapidly converging to 
the lowest possible variational energies. This ansatz is expected to provide a substantial part of the 
correlation energy for general complex atomic and molecular systems. The antisymmetric geminal 
term can be written as a single determinant even in the polarized cases. In general, therefore, the 
computational effort to sample this correlated wavefunction is not very demanding, the scaling of the 
algorithm with the number of atoms being comparable with the simplest Hartree-Fock calculation. 

We applied this Jastrow-geminal approach to atoms up to Z — 15, always getting good variational 
energies, by particularly improving those with a strong multiconfigurational nature. Our wavefunc- 
tion is very useful for Monte Carlo techniques, such as Fixed node. Indeed, the nodal surface 
obtained within this approach can be substantially improved through the geminal expansion. 

PACS numbers: 31.10,+z, 02.70.Ss, 02.70.Uu, 31.15.Ar, 31.15.Pf, 31.25.Eb, 71.10.Li 



I. INTRODUCTION 



One of the main goals in electronic structure calcu- 
lations is to deal with a wavefunction both accurate to 
predict the physical properties of a quantum system and 
simple enough to allow feasible computations of them. In 
particular, although multideterminantal CI like methods 
could be in principle very accurate, in practice the de- 
terminantal expansion becomes heavier and heavier from 
the computational point of view, as the number of de- 
terminants dramatically increases with the complexity of 
the electron system. On the other hand, a single determi- 
nantal wavefunction, kernel of methods like Hartree-Fock 
(HF) and Density Functional Theory (DFT), is some- 
times not sufficient to describe strongly correlated sys- 
tems, as for instance the transition metal compounds and 
the near degenerate shell structure of some atoms. 

Since 50's, the intensive efforts to explain theoretically 
the superconductivity have been highlighting the role of 
pairing in the electronic structure. The BCS wavefunc- 
tion belongs to an original ansatz in which the correlation 
is introduced through the product of pairing functions (in 
this context called Cooper pairs), already exploited in 
quantum chemistry by the pioneering work of Hurley et 
al. [l| to treat correlation effects in molecular properties. 
Their wavefunction was called antisymmetrized geminal 
power (AGP) that has been shown to be the particle- 
conserving version of the BCS ansatz 0|. It includes the 
single determinantal wavefunction, i.e. the uncorrelated 
state, as a special case and introduces correlation effects 
in a straightforward way, through the expansion of the 
pairing function (geminal): therefore it was studied as a 
possible alternative to the other multideterminantal ap- 
proaches. Although this ansatz seemed so appealing, it 



led to some expensive optimization procedures 001 with 
numerical problems 5, 6] in particular when applied to 
large systems, and so it turned out to be non competitive 
with respect to HF and CI. 

We show in this paper that the use of Monte Carlo 
methods can overcome the previous difficulties in opti- 
mizing the AGP wavefunctions. Two of the most ap- 
pealing features of these techniques are the possibility to 
tackle in a smart way the many-body interacting prob- 
lem, having the freedom in the choice of the functional 
form of the wavefunction, and to implement very effi- 
cient projection algorithms, like diffusion Monte Carlo 
(DMC) and Green's function Monte Carlo (GFMC). The 
trial wavefunction used in these methods is obtained by 
multiplying an antisymmetric term (usually called Slater 
term) to a Jastrow factor, which correlates the electrons 
and takes into account the interelectron cusp conditions 
the true wavefunction must fulfill. The Slater term can 
be either HF or CI or AGP like. As already pointed out 
by Umrigar .7] , the rate of convergence of CI expansion 
is increased by the Jastrow factor, just because it allows 
the wavefunction to have the correct cusps, otherwise 
present only asymptotically in the linear combination of 
determinants. 

In this paper we study the AGP-Jastrow (AGP + J) 
wavefunction, applied to the atomic problem. The wave- 
function we consider is actually a Resonating Valence 
Bond (RVB) state, first investigated in 1973 by Ander- 
son 0, H to study the high-temperature superconduc- 
tivity and later applied by Bouchaud and Lhuillier [Io| 
in Monte Carlo calculations of liquid 3 He. It is remark- 
able that this kind of wavefunction can describe well the 
ground state of all the atoms studied here, in particu- 
lar those light atoms affected by nondynamic correlation. 
Moreover, the Jastrow part is crucial not only in reducing 
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the pairing expansion but also in inducing a significant 
improvement of the nodal surfaces of the wavefunction, 
once we optimize both the Jastrow and the determinan- 
tal part at the same time. The minimization of the total 
energy can be efficiently performed by using an optimiza- 
tion procedure based upon the Stochastic Reconfigura- 
tion (SR) method [nj, which allows us to determine a 
large number of variational parameters, both for the Jas- 
trow and the determinantal part. Due to the simplicity 
of the AGP + J wavefunction and its capability to take 
into account resonating Slater configurations, this wave- 
function is expected to be appropriate also for molecules 
and more complex systems. 

The paper is organized as follows: in Sec. [H]we de- 
fine the wavefunction, in Sec. IIIII we describe the energy 
minimization method and finally in Sec. IIVI and we 
present detailed results and conclusions. 



We are going to extend the definition of the geminal 
wavefunction to a polarized system, i.e. a system with 
a different number of electrons for each spin. This gen- 
eralization of the geminal model was first proposed by 
Coleman [l5| and called GAGP. Without loss of gen- 
erality, one can assume that the number of up particles 
(N') is greater than the down ones (N*). In order to 
write down the many-body wavefunction with the gemi- 
nals, one needs to introduce — single particle spin 
orbitals $j not associated with any geminal and holding 
unpaired electrons. Once again one recovers the compact 
notation 10} for the spatial part of $agp (see Appendix) , 
but this time Aij is a x matrix defined in the fol- 
lowing way: 



A^ — 



,(r!) 



for j = l,N l 

for j = N L + l,JVf 



(6) 



II. FUNCTIONAL FORM OF THE WAVE 
FUNCTION 

The wavefunction we used in our QMC calculations is 
the antisymmetrized product of geminals (AGP), multi- 
plied by a simple symmetric Jastrow factor: 



where the index i ranges from 1 to 7V T . 

The pairing function can be expanded over a basis of 
single particle orbitals fig : 



M 

E 

i=l 



A^(r^*(r-L), 



(7) 



*(ri, ...,r N ) = # agp (ri, ■ ■ ■,r N )J{r 1 , . . .,r N ), (1) 

where the origin of the reference frame is the nuclear 
position. The two parts of the wavefunction Q will be 
described in detail below. 



A. The determinantal part 

Let $ be the pairing function (geminal) which takes 
into account the correlation between two electrons with 
opposite spin. If the system is unpolarizcd and the state 
is a spin singlet, the AGP wavefunction is 

*A G p(ri, ...,r N ) = A[9(t\, r|)$(r|, rj) • • • ^ N _ v r l N )}, 

where A is an operator that antisymmetrizes the product 
in the square brackets and the geminal is a singlet: 

< j, (r T ;r l ) = 0(r T jr i ) _L ( | U) _ UT))i (3) 

implying that 0(r,r') is symmetric under a permutation 
of its variables. Given this conditions, one can prove fhH 
that the spatial part of the ^agp can be written in a 
very compact form: 



*AGp(ri, ■• • ,rjv) = det(A i:j ), 



where A^ is a ^ x y matrix defined as: 



^•=<KrJ,rj). 



(4) 



(5) 



where ipi are general real or complex normalized func- 
tions and Xi are variational parameters that may depend 
on the chosen spatial symmetry of the geminal. Here- 
after, for simplicity, we do not assume that the orbitals 
are mutually orthogonal. 

For the application to atoms, in order the wavefunction 
has a definite angular momentum, it is convenient that 
the geminal is rotationally invariant around the nucleus. 
This requirement is achieved by taking the generic orbital 
ipi to be an eigenfunction of the single particle angular 
momentum operators I 2 and l z ; hence, the orbital will be 
denoted by: 



Ipnlm (r, 9,<t>) = R n i (r)Yi m [9, i 



(8) 



where Yi m are spherical harmonics with standard nota- 
tions and the radial part R n i depends on the principal 
quantum number n. In this way, the atomic geminal 
function takes the form: 



<XrV)=]T A^-irVw(r T )C im (r X ). (9) 

nl m— — l 

In the polarized case, the remaining orbitals ipj may 
change the total angular momentum and spin quantum 
numbers with the same rules valid for Slater-type wave- 
functions. Within our ansatz it is therefore possible to 
have definite total spin and angular momenta at least in 
all cases when the conventional Slater function does. 

The minimal geminal expansion is for M = N l , then 
the AGP wavefunction is reduced to the HF one; instead, 
if M is greater than N+, one can prove that the AGP is 
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equivalent to a linear combination of many Slater de- 
terminants. Therefore, within this functional form, one 
is able to introduce the correlation in the determinantal 
part of the wavefunction just by adding some terms in 
the geminal expansion. For a two-electrons closed-shell 
system, like Helium atom, or an ensemble of such non- 
interacting systems, this wavefunction is exact. For the 
other atoms, we will show that the AGP + J wavefunc- 
tion can lead to very good variational energies, even in 
the cases where the HF approximation is rather poor, 
especially for light elements. 

In order to optimize the radial part R n [ of the the single 
particle orbitals, we expand these radial functions in a 
Slater basis, in close analogy with Roothaan-Hartree- 
Fock calculations ^3]> namely using functions of the type: 



with n > 1, taking in principle as many different z^'s as 
required for convergence. 

In the Roothaan-Hartree-Fock the coefficients of the 
linear combinations are more involved, since the orthog- 
onality among all single particle states is required. In 
the Monte Carlo approach we have found that it is much 
simpler and more efficient to deal with non-orthogonal or- 
bitals, without spoiling the accuracy of the calculation. 
In fact, for light elements with Z < 15, studied here, it 
is possible to obtain almost converged results by using 
only two exponentials for each radial component (double 
zeta) . 

Hence, our single particle orbitals read in general 



R nl {r) =CV"- 1 (e- Zir + 



pe 



(11) 



where p is another variational parameter and C is the 
(10) normalization factor for the radial part R n i\ 



C 



//r, \-(2n+l) . o / , \-(2n+l) , 2/r, \-(2ra+l)Y 
l(2zi) >+2p{z 1 + z 2 ) ' +p (2z 2 ) v ' J 



(12) 



Actually p is not free for all the orbitals: indeed, for 
a more accurate variational wavefunction, we impose the 
electron-nucleus cusp condition 0,0, which is satisfied 
by the exact many body ground state and implies that 
each orbital must fulfill the following relation: 



dtj' 

dr 



(13) 



at r = (the hat denotes the spherical average). That 
condition is automatically obeyed by all but Is and 2s or- 
bitals of the type given in Eci llll Instead, the parameter 
p of Is orbital must be: 



P = 



Z\ — Z 

Z2 



(14) 



and for the 2s state, in order to fulfill Eq. (|13f) . we choose 
a functional form of the type: 



+ (p + ar)e 



(15) 



where a is a further variational parameter and p is given 
by 



P = 



Z 



Z2 



(16) 



In our study, we found that the presence of the a term 
leads to a very marginal improvement of the variational 
wavefunction, therefore we set a = and we kept the 
Is and 2s orbitals to have the same functional form, in 
order to reduce at most the total number of parameters. 



B. The Jastrow factor 

The Jastrow factor in our wave function is very sim- 
ple and has been widely used in previous Monte Carlo 
electronic structure calculations: 



J(ri, . . . ,r N ) = J|cxp(/(ri,rj)) 



(17) 



where the product is over all pairs of electrons and for 
simplicity the function / depends only on their relative 
distance r*y- and their spins 0$ and o~j, namely: 



ftwo body\Tij) 



(18) 



The value of a CTi(Tj . is fixed by the cusp condition at the 
coalescence points of two electrons and b aitTj contains at 
most three free variational parameters, as 6fj = is im- 
plied by the spatial symmetry of the Jastrow factor. The 
cusp conditions for parallel and antiparallcl spin electrons 
are different and yield two different values of a CTiCT . : 



1 



for o-j 



a = — for Gi 7^ <7j . 



(19a) 



(19b) 



As pointed out in Ref. 20, whenever a aiaj or b aiaj de- 
pend on the electron spins Oi ov, the wavefunction will 
be spin contaminated, i.e. it will not be an eigenstate of 
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the total spin operator S 2 . We have chosen to preserve 
the correct spin symmetry of the total wavefunction, by 
keeping a CTjCT . =1/2 and b aitTj — 6, hence fulfilling only 
the first condition (|19afl . Indeed, the cusp condition for 
electrons with parallel spins is much less important be- 
cause their probability to get close is clearly small, due 
to the Pauli principle. 

For the nitrogen atom, we checked the quality of this 
wavefunction with respect to a spin contaminated one 
with two variational parameters, frff = = b\ and 
^Ti = ^2, both of them reported in Table llTTl listed with 
the other atoms. In both the cases, we kept the gemi- 
nal expansion to be minimal (H F + J like wavefunction) . 
As reported in Table [I] the improvement in energy ob- 
tained by contaminating the spin wavefunction is rather 
negligible, and disappears when the FN DMC simulation 
is carried out. This implies that it is possible to obtain 
almost optimal nodes, without spoiling the spin symme- 
try and by using only one variational parameter for the 
Jastrow factor. 

In order to study the effectiveness of the Jastrow for 
a more accurate determination of the nodal surface, we 
have implemented a more involved Jastrow factor, includ- 
ing three-body correlation terms, and we have applied it 
to few atoms (Be and Mg). Indeed the Jastrow term / 
in Eq. El is a general function of the positions of two 
electrons and it has been parametrized similarly to the 
geminal function l(7)l. but truncated up to the I = 1 an- 
gular momentum: 



fir. 



It L J J 



ftwo bodyvij) 

i> (ri)Mrj) +?i(r<) (20) 



The functional form for the s-wave (ipo) and p-wave 
components can be chosen among different types; the 
most widely used in the literature are the polynomial 
[2lL I23 and the gaussian-polynomial mixed form [^J. 
In this work, we have selected the expansion over a gaus- 
sian basis, as reported in Table IVl This parametrization 
of the Jastrow factor, though being certainly less general 
compared with the best possible ones [23, includes the 
most significant part of the three-body correlation pflj . 
which involves two electrons and the nuclei. Our pur- 
pose, in fact, is to check whether it is possible to lower 
significantly the energy of the AGP + J wavefunction, 
whenever the Jastrow part of the wavefunction is system- 
atically improved together with the determinantal part. 



III. METHOD 



A. Minimization method 



of p variational parameters {a < l}k=i,...,p- Consider now 
a small variation of the parameters a% = a® + Sak . The 
corresponding wavefunction is equal, within the 

validity of the linear expansion, to the following one: 



(21) 



fe=i 



Therefore, by introducing local operators defined on each 
configuration x = {1*1, . . . , rjy} as the logarithmic deriva- 
tives with respect to the variational parameters: 



O k (x) 



d 

da k 



ln* T (x) 



and for convenience the identity operator O 
write ty'rp in a more compact form: 



(22) 



1, we can 



E 

k=0 



Sa k O k \^ T ) 



(23) 



where \^t) — \^t(o!q)) and Sao = 1- In general, Sao 
1, in that case the variation of the parameters will be 
scaled 



Sa k 



Sa k 
8a 



(24) 



and ty'j, will be proportional to ^^(a) for small j^. 

Our purpose is to set up an iterative scheme to reach 
the minimum possible energy for the parameters a, ex- 
ploiting the linear approximation for '^(a), which will 
become more and more accurate close to the conver- 
gence, when the variation of the parameters is smaller 
and smaller. We follow the stochastic reconfiguration 
method and define 



\*' t )=Psr(A-H)\* t ) 



(25) 



where A is a suitable larg e shift, allowing *£' T to have a 
lower energy than \Pr [H|> and Psr is a projection oper- 
ator over the (p + l)-dimensional subspace, spanned by 
the basis {Ok\^T)}k=o,...,p, over which the function \^' T ) 
has been expanded (Eq. 123(1 . Indeed, in order to deter- 
mine the coefficients {Sak}k=i,...,p corresponding to 
defined in Eq[2U one needs to solve the SR conditions: 

(* T |O fc (A - H)YS> T ) = (^ T \O k \%) for k = 0, . . . ,p 

(26) 

that can be rewritten in a linear system: 



1 s 



Ik 



(27) 



We have performed the wavefunction optimization by 
using the stochastic minimization of the total energy 
based upon the Stochastic Reconfiguration ( SR) tech- 
nique, already exploited for lattice systems [Tj. Let 
^/^(a ) be the wavefunction depending on an initial set 



where ""' 



is the covariance matrix and 

Ik 



s- = (m T \O l O k \y T , 
f k = {^ T \O k (A - H)\^ T ) is the known term; both s 
and f k are computed stochastically by a Monte Carlo 
integration. These linear equations (|27jl are very similar 
to the ones introduced by Filippi and Fahy 26] for the 
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energy minimization of the Slater part. In our formula- 
tion, there is no difficulty to optimize the Jastrow and 
the Slater part of the wavefunction at the same time. 

After the system l|27|l has been solved, we update the 
variational parameters 

ak=a f + pL for k=l,...,p (28) 
Sa 

and we obtain a new trial wavefunction ^>t(ck). By 
repeating this iteration scheme several times, one ap- 
proaches the convergence when -> for k ^ 0, and 
in this limit the SR conditions l|26(l implies the Eulcr 
equations of the minimum energy. Obviously, the solu- 
tion of the linear system (|27J) is affected by statistical 
errors, yielding statistical fluctuations of the final vari- 
ational parameters even when convergence has been 
reached, namely when the {ctk}k=l ...,p fluctuate without 
drift around an average value. We perform several itera- 
tions in that regime; in this way, the variational param- 
eters can be determined more accurately by averaging 
them over all these iterations and by evaluating also the 
corresponding statistical error bars. 

It is worth noting that the solution of the linear system 
(|27|l depends on A only through the Sao variable. There- 
fore the constant A indirectly controls the rate of change 
in the parameters at each step, i.e. the speed of the al- 
gorithm for convergence and the stability at equilibrium: 
a too small value will produce uncontrolled fluctuations 
for the variational parameters, a too large one will al- 
low convergence in an exceedingly large number of iter- 
ations. The choice of A can be controlled by evaluating 
the change of the wavefunction at each step as: 

'\~* T ' 2 = E ^ Sa » (») 

' T ' k,k'>0 

By keeping this value small enough during the opti- 
mization procedure, one can easily obtain a steady and 
stable convergence. 

Finally, we mention that the stochastic procedure is 
able in principle to perform a global optimization, as dis- 
cussed in Ref. [n] for the SR and in Ref. for the 
Stochastic Gradient Approximation (SGA), because the 
noise in the sampling can avoid the dynamics of the pa- 
rameters to get stuck into local minima. 



B. Variational and Diffusion Monte Carlo 

We performed standard Variational (VMC) and Diffu- 
sion Monte Carlo (DMC), the latter within the so called 
Fixed Node (FN) approximation, which allows to obtain 
the lowest energy state with the same nodes of a trial 
wavefunction. As trial state for FN, we have used the 
VMC wavefunction optimized with the SR method de- 
scribed in the previous section. 



IV. RESULTS 

We have carried out Quantum Monte Carlo calcula- 
tions for atoms from Li to P, using the antisymmetrized 
geminal power times the Jastrow factor (AGP + J) to 
describe the atomic electronic structure. We performed 
the optimization of both the geminal and the Jastrow 
part minimizing the energy with the method described 
in Sec lIIII For all the atoms, we considered the minimal 
geminal expansion, corresponding to the HF single de- 
terminant, and the simplest Jastrow factor with a single 
parameter (|18|) , reported in Table 11111 To improve the 
antisymmetric part, we increased the number of orbitals 
in the geminal expansion, and for Be and Mg atoms 
we also systematically considered an improved Jastrow 
term, such as the three-body one described above (see 
Tables llVl and IV |) . As one can notice from the Tables, 
our wavefunction parametrization is very compact, even 
in the case of the mostly correlated states, since it con- 
tains always a relatively small number of parameters for 
each atom. 

In order to judge the outcome of our calculations, we 
computed the correlation energies and in particular its 
fraction with respect to the exact ground state energy 
for the nonrelativistic infinite nuclear mass Hamiltonian, 
estimated in Ref. The quality of the variational wave- 
function can be seen by computing the expectation value 
of the energy by means of the VMC calculations. Fur- 
thermore, we carried out DMC simulations within the FN 
approximation, which allows to optimize the amplitude of 
the wavefunction inside each nodal volume, where its sign 
is given and fixed by the variational state. Therefore, the 
DMC energy depends on the quality of the nodal struc- 
ture of the variational wavefunction and the capability 
of improving the nodes during the optimization is cru- 
cial to obtain almost exact DMC energy values. To that 
purpose, it is very important to have a variational func- 
tional form appropriate to reproduce the correct nodes. 
We show that the AGP + J wavefunction satisfies well 
this requirement, yielding in all the atoms studied here 
very good DMC results. The VMC and DMC energies 
are listed in Table U in Figures ^ an d El we plot the per- 
centage of the correlation energy recovered respectively 
by VMC and DMC calculations for different atoms and 
wavefunctions. 

The VMC calculations with the minimal geminal ex- 
pansion and the two body Jastrow factor yield from 60% 
to 68% of the total correlation energy, with the exception 
of the Li atom, where 91.4% of the correlation energy 
is obtained. Therefore, there is a sizable loss of accu- 
racy in going from Li to Be, the worst case being the 
Boron atom. The corresponding DMC simulations get a 
large amount of the energy missing in the VMC calcula- 
tions, recovering from 87.7% to 99.9% of the total cor- 
relation energy, but the dependence on the atomic num- 
ber shows the same behavior found in VMC: the worst 
results are obtained for Be, B and C atoms, due to the 
strong multiconfigurational nature of their ground states. 
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As well known, one can improve substantially the varia- 
tional state of those atoms including not only the Is 2 2s 2 
configuration but also the 2s 2 2p 2 , because of the near de- 
generacy of 2s and 2p orbitals. In this case the AGP + J 
ansatz is particularly efficient: by adding just one term 
in the geminal expansion, we are able to remove this loss 
of accuracy in the correlation energy both in the VMC 
and the DMC calculations. 

In Table ^ we summarize some results obtained for 
Be in previous works and compare them with ours. AGP 
calculations of atoms have been performed only few times 
so far, the best one for Be is reported in the last row of 
the Table HU Kurtz et al. were able to recover 84% 
of correlation energy using a geminal expansion with a 
very large basis; our variational AGP + J wavefunction 
reaches 94% with a much smaller basis (Is, 2s and 2p or- 
bitals). By including a three-body Jastrow factor, 97.5% 
of the correlation energy is finally obtained, which is com- 
parable with the best multiconfigurational wavefunctions 
previously studied [24|. 

This outcome highlights the importance of the Jastrow 
in reducing the geminal expansion and yielding a better 
energy. The nodal surface can be substantially improved 
with the present approach, because the pairing expan- 
sion contains implicitly not only the four determinants 
ls 2 2s 2 and ls 2 2p 2 , but also the other three 2s 2 2p 2 and 
six 2p 2 2p 2 , which can improve further the wavefunction. 
Indeed, the geminal expansion reduces exactly to four 
determinants in the limit \2 S — > and A2 P — > with 
constant ratio ^ 2jL . The fact that the minimum energy is 

obtained for A2 S ^ and \2 P ^ (see Tables Hvl and Ivjl 
clearly shows that the energy can be lowered by consid- 
ering the remaining configurations described above. In- 
deed our DMC energies are slightly lower than the ones 
by Huang et «/.|2~H|. to our knowledge the best available 
ones obtained with the four determinants: one ls 2 2s 2 
and three ls 2 2p 2 . In order to determine accurate nodes 
for the corresponding DMC calculation they used the 
two-body Jastrow factor similar to the one H18|) we used 
or an highly involved three-body term much more com- 
plete than ours (for this reason our corresponding VMC 
energy is slightly worse in this case). We also verified, 
therefore, that a more accurate description of the Jas- 
trow factor (which do not affect directly the nodes) is 
crucial to obtain better nodes, provided the variational 
parameters, belonging to both the Jastrow and the de- 
terminantal part, are optimized altogether. For instance, 
in Ref. |2^ the authors optimized only the coefficients in 
front of the four determinants ls2s — ls2p and not the 
orbitals, obtaining for Be energy not comparable with 
the best possible ones. The AGP + J is simple enough to 
allow a feasible parametrization of the variational state, 
by capturing the most important correlation. 

We found that also Mg, Al and Si, the second row 
atoms corresponding to Be, B and C in the first row, 
have a quite strong multiconfigurational character, in- 
volving here 3s and 3p orbitals. Analogously to the 
Beryllium case, for the Mg we have optimized both the 



two-body and the three-body Jastrow factor, together 
with the AGP wavefunction containing 3s-3p resonance. 
In this case, although at the variational level the three- 
body wavefunction is much better than the two-body one 
(see Fig. ^| , that difference disappears almost completely 
in the DMC results. This shows that the correction of 
the nodal surface allowed by the more accurate three- 
body Jastrow does not seem to be crucial as in the Be 
atom. On the other hand, the effect of the AGP expan- 
sion is significant in improving further the DMC calcula- 
tion, which already yields good FN energies even with a 
simple HF + J trial wavefunction for atoms heavier than 
C (percentage of correlation energy always greater than 
92%). By adding the 3p contribution to the geminal we 
were able to recover 96.8% of the correlation energy of 
Mg (see Fig. Also for Al the presence of the 3p or- 
bital is significant in reducing the DMC energy, and for 
Si it seems important only in the VMC value. 

Finally, by using the AGP + J wavefunction, we opti- 
mized some atoms (from N to Na) not affected by non- 
dynamic correlation; here, in order to obtain an improve- 
ment in the VMC and in the DMC energies, we needed 
a bigger basis (3s2pld) to be used in the geminal expan- 
sion. 



V. CONCLUSIONS AND PERSPECTIVES 

In this work we have introduced a variational wave- 
function which contains the main ingredients of electron 
correlations: the Jastrow factor, that allows to satisfy the 
electron-electron cusp condition, and the geminal expan- 
sion, that allows to consider a correlated multiconfigura- 
tion state, with a numerically feasible scheme, namely by 
evaluating only a single but appropriately defined deter- 
minant. 

The application to atoms is particularly successful for 
low atomic number, where Hartree-Fock is particularly 
poor, due to the almost degenerate 2s — 2p shells. The 
case of Beryllium is an interesting benchmark. Indeed, 
by considering the change of the geminal part altogether 
with the Jastrow term, we obtained an excellent repre- 
sentation of this correlated atom. Our results, presented 
in Table [H] are not only comparable but appear even 
better than the best multideterminantal schemes (using 
e.g. four Slater determinants), showing that it is possible 
to represent non trivial correlated states by properly tak- 
ing into account the interplay of the Jastrow term and 
the determinantal part of the wavefunction. Our vari- 
ational energies for the other atoms (see Table can 
be substantially lowered because we have considered, in 
this first application, a wavefunction with the two-body 
Jastrow factor. 

As well known the variational energy of the Hartree- 
Fock wavefunction cannot be improved by extending 
the variational calculation to a larger basis including 
all particle-hole excitations applied to the Hartree-Fock 
state. Analogously, the geminal wavefunction is not only 
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stable with respect to these particle-hole configurations 
, but also to all possible states obtained by destroy- 
ing a singlet pair on some orbital and creating another 
one on another orbital. Though this wavefunction can 
take into account a big number of configurations which 
may allow an energy improvement, obviously it cannot 
include everything within a single geminal, otherwise 
the wavefunction would be exact. Indeed there exist 
multiconfigurational states that are known to be impor- 
tant for atoms like C or those with larger Z [29|. and 
that involve complicated multi-particle excitations to the 
Hartree-Fock state. These ones cannot be reduced to cre- 
ation/destruction of singlet pairs and therefore cannot be 
handled with a single geminal function. However in our 
study we have found that the single geminal function 
with the proper Jastrow factor already provides satisfac- 
tory results for all atoms, yielding more than 93% of the 
correlation energy in all cases studied by carrying out 
DMC simulations. 



The extension of this approach to molecules or more 
complex electronic systems is straightforward, and is in- 
deed particularly promising. As pointed out in Ref. |2?| 
the geminal wavefunction for a diatomic system can cor- 
rectly describe the interatomic Born-Oppenheimer po- 
tential from small to large distances, where, in this limit, 
the wavefunction of isolated atoms can be smoothly re- 
covered. This important property cannot be satisfied 
within the Hartree-Fock theory, even for the simplest H2 
molecule (without contaminating the singlet ground state 
wavefunction) . 



For an electronic system with many atoms, the gem- 
inal expansion together with the Jastrow term is very 
similar to the so called Resonating Valence Bond (RVB) 
expansion, introduced by Pauling (see e.g. 30]) long time 
ago, and revived recently by P.W. Anderson to consider 
the properties of strongly correlated electronic systems. 
The geminal part, when expanded in terms of Slater de- 
terminants, yields a very large and non trivial number 
of configurations, which increases exponentially with the 
number of atoms considered. The Jastrow factor in this 
case suppresses the weight of those configurations with 
two electrons close to the same atomic orbital, correctly 
describing the effect of the strong Coulomb repulsion. We 
see therefore the remarkable advantage of this approach. 
Just for complex systems with many atoms an RVB wave- 
function corresponding to an exponentially large number 
of configurations can be efficiently used for a more ac- 
curate description of electron correlation. Contrary to 
the conventional RVB expansion, it is appealing , not 
only from the computational point of view, that these 
properties can be obtained by sampling a single deter- 
minant wavefunction within the Quantum Monte Carlo 
techniques. 
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APPENDIX: SPIN POLARIZED GEMINAL 
WAVEFUNCTION 

In this appendix, we consider the most general gem- 
inal wavefunction with definite spin S = — 2 — , where 
2Vf (N± ) is the number of spin-up (spin-down) electrons 
and Nf > is assumed. To this purpose we introduce 
second quantized fermionic fields (see e.g. Fetter and 
Walecka |U) ip^(r,a) and ip(r,a), where r is the elec- 
tron position and a = ±1/2 is its spin projection along 
the z-axis. These fields satisfy the canonical anticommu- 
tation rules: 

{^(r, a), ^t(r', a')} = 5 aa ,S(r - r'). (A.l) 

In these notations, the most general wavefunction with 
definite spin can be formally written in the following way: 

I*} =Pn 1] 4a exp(*t)j ), (A.2) 

i=N l + l 

where Pjv is the projection on the given number of par- 
ticles N = + Ni , |0) denotes the vacuum of electrons 
and tpj 1 is the most generic (Bogoliubov) orbital with 
spin S = 1/2: 

i'U = J d * (<^(r)VM) +^ > (r)V t (r,T)) , (A.3) 

which is defined by the orbital functions 4>f for the cre- 
ation of a particle with spin up and ipf for the annihi- 
lation of a particle with spin down. For instance, a con- 
ventional Slater determinant of spin-up particles can be 
written as Yl i ^j||0), where <frf — 0. It is clear therefore 
that this representation is more general and may provide 
a wavefunction ^ much more reach than the conventional 
Slater determinants. 

Finally, the pairing creation operator $t j s a singlet, 
namely exp($t)|o) has spin zero, and is defined by a 
generic symmetric function $(r, r') = $(r',r): 

$t = J dv J dv' $(r',r)V' t (r,|)V' t (r',T). (A.4) 

Our purpose is to show here that the value of the 
wavefunction W can be simply computed, similarly to 
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a conventional Slater determinant, on each configuration Indeed the value F of the wavefunction on (a;| is: 
x = {r^f, . . . , rjVi ,J.}> where r^j are the positions of spin- 
up particles and are the spin-down ones. These con- 
figurations can be generally written as: 

(xi=<oin^,T)n^->i)- ( a - 5 ) 

i=l j=l 

I 



F = (x|#) = ( 



Ni 



(A.6) 



Now we insert the identity exp(— ex p($t) between 
each fermionic field in the above equation l|A.6(l : 



F = (0 |exp($ t )exp(-$ t )-0(ri,T)exp($ t ) 
• • • exp(-$ t )V'(rjv i , I) exp($ t ) • • • 



•exp(-$ t )Vj VTiT exp($ t ) 



(A.7) 



Exploiting the relation valid for generic operators A, B 
and C: 

exp(-A) B exp(A) = B-[A,B} + ^ [A, [A, B}]+... 

(A.8) 

one is able to evaluate the following terms: 

exp(-$ t )V'(r i ,T)cxp($ t ) = 
^(r is T)- Jdr $(r iit ,r)^( r ,|) 

exp(-$ t )V'(r l ,4)exp($ t ) = 
^(r,,4)+ Jdr <S>{r,r lA )^(r,^ 

cxp(-$ t )V'j t T exp($ t ) = 
4 A + J dr J dr ' *(r,r')</'<(r , )V' t (r,T) (A.9) 
In order to derive the above relations, notice that all the 



terms in the RHS of Eq. IA.8I are always zero beyond 
the first two. After substituting the expressions in Eq. 
lA~7land by using (0|exp($t) = (o|, V(r,cr)|0) = and 
(0\i/j'(r,a) = 0, one can iteratively apply the canonical 
commutation rules (|A.1|I and a simplified result for F is 
obtained: 



F = ( 



N, 



n^ r *>t)n^,i 



(A.10) 



where ^ is the creator of an orbital function of the type 
(|A.3fl . with transformed orbitals: 



i(r) 
i(r) 



$(r,r a ) (A.lla) 
fori = l,-- - ,Ni 

0>(r)+ Jdr' $(r,r')<^(r') (A.llb) 
for i = JV; + !,-•• ,JV T 



Then the final value of F can be simply computed by 
a single determinant, as it represents just the value of a 
iVf x JVj- Slat er determinant with orbitals given in HA. Hall 
and (|A.llb|l on the spin- up configurations, yielding the 
final expression © reported in the text. 
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TABLE I: Total energies in variational (Evmc) an d diffusion (Edmc) Monte Carlo calculations; the percentages of correlation 
energy recovered in VMC(EY MC (%)) and DMC (E° MC (%)) have been evaluated using the exact (E ) and Hartree-Fock 
(Ehf) energies from RefHa M is the number of terms in the pairing expansion. The energies are in Hartree. 





M 


E 


Ehf 


Evmc 


eY mc (%) 


Edmc 


E? MC (%) 


Li 


1 


-7.47806 


-7.432727 


-7.47415(10) 


91.4(2) 


-7.4780(2) 


99.9(4) 


Be 


2 


-14.66736 


-14.573023 


-14.63145(5) 


61.9(5) 


-14.6565(4) 


88.5(4) 




5 






-14.661695(10) 


93.995(11) 


-14.66711(3) 


99.73(3) 




5 a 






-14.66504(4) 


97.54(5) 


-14.66728(2) 


99.92(2) 


B 


2 


-24.65391 


-24.529061 


-24.6042(3) 


60.3(2) 


-24.63855(5) 


87.7(4) 




5 






-24.62801(4) 


79.25(4) 


-24.6493(3) 


96.3(3) 


C 


2 


-37.8450 


-37.688619 


-37.7848(6) 


61.5(4) 


-37.8296(8) 


90.2(5) 




5 






-37.7985(7) 


70.3(4) 


-37.8359(8) 


94.2(5) 


N 


2 


-54.5892 


-54.400934 


-54.52180(15) 


64.20(8) 


-54.57555(5) 


92.7(3) 




2'' 






-54.52565(15) 


66.20(8) 


-54.5753(4) 


92.6(2) 




14 






-54.5263(2) 


66.62(11) 


-54.5769(2) 


93.47(10) 





3 


-75.0673 


-74.809398 


-74.9750(7) 


64.2(3) 


-75.0477(8) 


92.4(3) 




14 






-74.9838(6) 


67.6(2) 


-75.0516(9) 


93.9(3) 


F 


4 


-99.7339 


-99.409349 


-99.6190(8) 


64.6(3) 


-99.7145(15) 


94.0(5) 




14 






-99.6315(7) 


68.4(2) 


-99.7141(6) 


93.91(18) 


Ne 


5 


-128.9376 


-128.547098 


-128.8070(10) 


66.6(3) 


-128.9204(8) 


95.6(2) 




14 






-128.8159(6) 


68.83(17) 


-128.9199(7) 


95.47(18) 


Na 


5 


-162.2546 


-161.858912 


-162.1334(7) 


69.37(19) 


-162.2325(10) 


94.4(2) 




o 
y 






1 fiO 1/1^/1 (1\ 


71 Q1 (A «\ 


i o^vfVi r\\ 

— IDZ.Zo / U^lUj 


yo.o^z ) 


Mg 


6 


-200.053 


-199.614636 


-199.9113(8) 


67.67(19) 


-200.0327(9) 


95.4(2) 




9 






-199.9363(8) 


73.38(19) 


-200.0375(10) 


96.5(2) 




9" 






-200.0002(5) 


87.95(12) 


-200.0389(5) 


96.79(11) 


Al 


6 


-242.346 


-241.876707 


-242.1900(9) 


66.77(19) 


-242.3215(10) 


94.8(2) 




9 






-242.2124(9) 


71.53(19) 


-242.3265(10) 


95.8(2) 


Si 


6 


-289.359 


-288.854363 


-283.1875(10) 


66.0(2) 


-289.3275(10) 


93.8(2) 




9 






-289.1970(10) 


67.9(2) 


-289.3285(10) 


94.0(2) 


F 


6 


-341.259 


-340.718781 


-341.0700(10) 


65.0(2) 


-341.2220(15) 


93.2(3) 



"Wavefunction with a three body Jastrow factor. 
'Wavefunction with a two body Jastrow factor spin contaminated. 



TABLE II: Comparison of the energies obtained by various authors for Be. 



basis Jastrow VMC DMC 

present work 2sTp two body -14.661695(10) -14.66711(3) 

Huang et al. [24] 2slp two body -14.66088(5) -14.66689(4) 

present work 2slp three body -14.66504(4) -14.66728(2) 

Huang ETA.\2£ 2slp three body -14.66662(1) -14.66723(1) 

Sarsa et aZ.[22] 2slp three body -14.6523(1) 

Kurtz et a/.[2JJ 6s3p2d none -14.6171 
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TABLE III: HF+J (two body) wavefunctions 

Parameters of the Jastrow and the pairing function with the notation described in the text. means that the corresponding 
parameter has to be evaluated from the cusp condition in Eq. 1151 The line over the orbitals label refers to the unpaired ones. 
The values are given with the statistic error due to the stochastic approach in the minimization. 





b 


orbital 


Z x 


z 2 


P 


A 


Li 


0.731(3) 


Is 


3.556(2) 


2.3741(15) 


# 


1.0 






2s 


1.4289(12) 


0.5380(2) 


# 




Be 


0.773(2) 


Is 


4.569(6) 


3.289(5) 


# 


1.0 






2s 


2.602(3) 


0.7850(8) 


# 


1.0 


B 


0.877(2) 


Is 


5.569(5) 


4.195(4) 


# 


1.0 






2s 


3.527(2) 


1.0788(2) 


# 


1.0 






2p 


2.437(4) 


1.1001(5) 


0.2664(3) 




C 


0.990(2) 


Is 


6.533(6) 


5.075(8) 


# 


1.0 






2s 


4.475(3) 


1.3552(3) 


# 


1.0 






2p 


2.9835(12) 


1.3886(5) 


0.2374(4) 




N 


1.110(3) 


Is 


7.461(8) 


5.901(15) 


# 


1.0 






2s 


5.387(3) 


1.6323(1) 


# 


1.0 






2^ 


3.4908(9) 


1.6373(4) 


0.2036(3) 




N 


contaminated 


Is 


7.2738(17) 


5.522(7) 


# 


1.0 




6 n =0.940(3) 


2s 


5.254(2) 


1.6692(15) 


# 


1.0 




& TT =0.624(4) 


2p 


3.4544(4) 


1.66559(13) 


0.21593(9) 




O 


1.152(2) 


Is 


8.303(5) 


6.53(3) 


# 


1.0 






2s 


6.428(2) 


1.9362(5) 


# 


1.0 






2p 


3.9921(7) 


1.84932(15) 


0.1648(2) 


1.0 






2p 


3.9921(7) 


1.84932(15) 


0.1648(2) 




F 


1.226(9) 


Is 


9.171(3) 


6.82(2) 


# 


1.0 






2s 


7.380(2) 


2.2402(4) 


# 


1.0 






2p 


4.5025(2) 


2.0758(2) 


0.1451(10) 


1.0 






2p 


4.5025(2) 


2.0758(2) 


0.1451(10) 




Ne 


1.321(3) 


Is 


10.103(4) 


7.01(4) 


# 


1.0 






2s 


8.341(10) 


2.5319(5) 


# 


1.0 






2p 


5.0273(10) 


2.3155(4) 


0.1358(2) 


1.0 


Na 


1.514(7) 


Is 


11.102(2) 


7.58(5) 


# 


1.0 






2s 


8.823(2) 


2.8795(4) 


# 


1.0 






2p 


5.7178(9) 


2.8092(8) 


0.18027(18) 


1.0 






3s 


1.540(3) 


0.734(2) 


0.1265(2) 




Mg 


1.654(6) 


Is 


12.0855(15) 


7.96(5) 


# 


1.0 






2s 


9.349(4) 


3.2611(16) 


# 


1.0 






2p 


6.3795(10) 


3.2905(5) 


0.21569(17) 


1.0 






3s 


1.9288(10) 


1.02214(13) 


0.1449(2) 


1.0 


Al 


1.812(8) 


Is 


13.0791(17) 


8.43(7) 


# 


1.0 






2s 


9.846(8) 


3.646(2) 


# 


1.0 






2p 


7.0560(15) 


3.7873(7) 


0.2597(3) 


1.0 






3s 


2.24(8) 


1.2627(2) 


0.1445(2) 


1.0 






3p 


1.83(2) 


0.8866(2) 


0.1012(7) 




Si 




Is 




8 77 (s\ 


ft 


1 n 






2s 


10.40(5) 


4.0275(8) 


# 


1.0 






2p 


7.703(10) 


4.261(5) 


0.2944(12) 


1.0 






3s 


2.468(5) 


1.46(2) 


0.1241(2) 


1.0 






3p 


2.274(2) 


1.16(2) 


0.1227(3) 




P 


2.074(9) 


Is 


15.053(4) 


9.10(15) 


# 


1.0 






2s 


10.997(5) 


4.4269(5) 


# 


1.0 






2p 


8.346(2) 


4.7519(12) 


0.3273(17) 


1.0 






3s 


2.7386(5) 


1.62(2) 


0.1129(3) 


1.0 






3p 


2.5918(18) 


1.3479(7) 


0.1125(2) 
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TABLE IV: AGP+J (two 


body) wavefunctions 


The notations 


are the same as 


in Table UTTI 




b 


orbital 




^2 


P 




Be 


1.009(5) 


Is 




Q Q1 1 <A(1\ 
O.O L LOyO } 


M- 

w 


i n 

1 .u 






2s 


9 97561 (1 6 s ) 


78Q55('2 S ) 


41 
TT 


— 1 95(81 10~ 3 






2p 


o.ooy y^ ^ 




f) 8S7S('1 8"! 
u.oo l oyLoj 


Q 4Qf 1 i n _ 


B 


1.005(7) 


Is 


c; fi^no/R^ 


A 9041 ^7^ 


M- 

ft 


1 






2s 




1 07904^9^ 


M- 
W 


4 41 ^7^ 1 — 3 






2p 


Q 701 ('7 S \ 


1 4(S9QS('1' : i s \ 




fi Q9M 0"! 1 n~ 






2p 




1 07107^4^ 


981 D9('8^ 




C 


1.036(5) 


Is 




^ 0404^^ 
o.uiy^^o ) 


4J- 


i n 






2s 


A Af\QA.(*\\ 

'-t. L ±\JiJ L ±\0 } 




M- 
ft 








2p 


^t.zoy yo ) 


1 ^^17^9^ 










2p 


9 Q^7^R(R\ 
£.ijo i ooyo j 


1 ^fi9SS('9^ 


fl 9'}97fi('9^ 




N 


1.124(5) 


Is 




O.OUU^Z J 


4J- 


1 






2s 


c; ciO^I (19^ 




M- 
ft 


— 9 ^A^(Q\ i n - 






2p 


o. / Uo^y ^ 












3s 






u.u 


7 ^91 /^4^ 1 - 






3p 


9 1 Q1 /I 

z. ±y i^o ^ 




u.u 


o.QoQyl-O) 1U 






3d 


9 fiAQRffi^ 




u.u 


1 SSI ^4^ 1 n~ 






2p 


q 4SQ77/'1 1 ^ 




fl 9fl c i^7^(' 1 9") 




O 
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"In the case of dynamic 


correlation, we 


used single zeta orbitals 









beyond the HF ones. 
''For Mg, Al and Si, we optimized only the A's related to the 3s — 
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TABLE V: AGP+J (three body) wavefunctions 

The notations are the same as in Table ITTT1 The parametrization for the three-body Jastrow factor is also reported. 



Be 



geminal 


Zx 


z 2 


P 


A 




Is 


7.48(8) 


3.545(2) 


# 


1.0 




2s 


2.403(3) 


0.8181(3) 


# 


3.95(9) 10~ 3 




2p 


2.86(5) 


1.1064(9) 


1.245(9) 


-6.27(9) 10~ 4 




Jastrow 2 body 


b 




0.7935(19) 


Jastrow 3 body 


Zx 


z 2 


a 


ax 


a 2 




2.238(4) 


15.9(2) 


-0.3491(4) 


0.1122(4) 


-0.6067(12) 




0.19775(14) 


6.060(10) 


0.16535(9) 


1.907(2) 




Mg 


geminal 


Zx 


z 2 


P 


A 




Is 


12.0 


0.0 


# 


1.0 




2s 


10.784(2) 


3.1685(4) 
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2.40(18) 10 _1 




■2p 


6.6199(4) 


3.25296(19) 


0.19285(9) 


8.56(2) 10" 1 




3s 


1.9220(3) 


1.01249(15) 


0.08312(8) 


1.674(13) 10" 4 




3p 


2.599(12) 


1.2738(3) 


0.535(14) 


-2.152(17) 10" 5 




Jastrow 2 body 


b 




1.191(5) 


Jastrow 3 body 


Zx 


z 2 


a 


ax 


a 2 


Mr) c 


1.465(6) 




-0.3919(5) 


1.3867(9) 


-0.891(2) 




13.03(3) 


1.074(10) 


0.8789(2) 


0.29822(5) 





a ipp{r) = a (exp(-Zir 2 ) + ax exp(-Z 2 r 2 ) + a 2 ) 
b ipx( r ) = r ao(exp(-Zir 2 ) + ax cxp(-Z 2 r 2 )) 
c -0 o (r) = a ((l + axr 2 ) exp(-Zir 2 ) + a 2 ) 



Variational Monte Carlo 
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FIG. 1: VMC correlation energies for HF + J 2 (minimal geminal expansion with a two-body Jastrow factor), AGP + J 2 (best 
geminal expansion) and AGP + J3 (best geminal with a three-body Jastrow factor) 



Diffusion Monte Carlo 
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FIG. 2: DMC correlation energies obtained by HF + ,h, AGP + ,h and AGP + Js wavefunctions 



